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The stochastic motion in a nonhomogeneous medium with traps is studied and diffusion prop¬ 
erties of that system are discussed. The particle is subjected to a stochastic stimulation obeying 
a general Levy stable statistics and experiences long rests due to nonhomogeneously distributed 
traps. The memory is taken into account by subordination of that process to a random time; then 
the subordination equation is position-dependent. The problem is approximated by a decoupling of 
the medium structure and memory and exactly solved for a power-law position dependence of the 
memory. In the case of the Gaussian statistics, the density distribution and moments are derived: 
depending on geometry and memory parameters, the system may reveal both the subdiffusion and 
enhanced diffusion. The similar analysis is performed for the Levy flights where the finiteness of the 
variance follows from a variable noise intensity near a boundary. Two diffusion regimes are found: 
in the bulk and near the surface. The anomalous diffusion exponent as a function of the system 
parameters is derived. 


I. INTRODUCTION 

The diffusion is anomalous when the variance rises with time slower or faster than linearly and the density distri¬ 
bution differs from the normal distribution. This means that the central limit theorem is violated. One can expect 
that the theorem is not valid for transport processes if memory effects are present, e.g. when a medium contains 
traps. Traps hamper the transport and, as a consequence, a subdiffusion emerges - as well as a stretched-Gaussian 
asymptotics of the distribution [l|. We encounter such a situation for disordered media with impurities and defects. 
If the waiting time distribution is not prescribed to a given position and changes at each visit of the particle, we 
are dealing with an annealed disorder which corresponds to a renewal process. Then the continuous-time random 
walk (CTRW) is well suited to model the anomalous diffusion. If, on the other hand, the particle exercises a space 
structure which slowly evolves with time, the trapping time at a given site is the same for each visit of this site and 
a correlation between trapping times emerges (a quenched disorder) [^. The quenched trap model Si takes into 
account that the particle may remember the rest time at a given point. Then the time-dependence of the variance can 
be derived by averaging over disorder i or using a renormalisation group approach S; it takes a form ^ 
where 0 < /? < 1 characterises the memory. 

The central limit theorem may not apply also in the Markovian case since a nonhomogeneous structure of the 
environment makes subsequent random stimulations mutually dependent. Considering, in particular, the diffusion on 
fractals one must take into account the self-similar medium structure; it is described by the Fokker-Planck equation 
with a variable diffusion coefficient @ (for a non-Markovian generalisation see [3j0)- CTRW involves, in general, 
a coupled jump density distribution. When the waiting-time distribution is position-draendent, the Fokker-Planck 
equation, corresponding to the master equation, contains the variable diffusion coefficient [9| and, in the non-Markovian 
case, a variable order of the fractional derivative M- For such inhomogeneous problems, CTRW implies a stretched- 
Gaussian shape of the density distribution and predicts the anomalous diffusion. 

If, on the other hand, the jump length does not obey the normal distribution but is governed by a general Levy 
stable distribution (Levy flights) the variance does not exist; such long jumps are frequently observed in many areas 
of science m- They may be directly related to a specific topology of the medium and then the transport description 
requires a variable diffusion coefficient: the folded polymers are a well-known example [T^ . If a composite medium 
consists of many layers, the fractional equation is complicated and contains position-dependent both the diffusion 
coefficient and the order of the fractional derivative [Q. However, presence of the Levy flights does not need to 
imply infinite fluctuations because any physical system is finite and, if one introduces a truncation of the distribution 
tail, the diffusion properties are well-determined. It has been recently demonstrated that cracking of heterogeneous 
materials reveals a slowly falling power-law tail of the local velocity distribution of the crack front [l^ but, despite 
that, the authors were able to determine the variance due to the finiteness of the system. 

Systems with memory can be conveniently handled in terms of a Langevin equation by introducing an auxiliary 
operational time. Process given by this equation is subsequently subordinated to the random physical time by means 
of a one-sided probability distribution with a long tail [l5|, [l6[ ; such a system of the Langevin equations is used as an 
alternative formulation of CTRW [ 1 I 0 and also applied to quenched random media [l^ where the intensity of the 
random time distribution depends on the position. CTRW was applied to transport in the heterogeneous media by 
using fractional derivatives and a subordinated system of Langevin equations . A subordination formalism, 
which directly takes into account that the memory in nonhomogeneous systems must depend on the position, has 
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been recently proposed . This dependence was introduced via a variable intensity of the random time distribution 
and models the influence of the trap geometry on the local time lag. The system is then described by a set of two 
Langevin equations, 

dx{T) = ri^dr) 

dt{T) = g{x)^{dT), (1) 


where a non-negative random time intensity g{x) defines a position-dependence of the memory effects and results from 
the presence of traps. Increments of the white noise g are determined, in general, by the a-stable Levy distribution, 
La{x), and ^(dr) stands for a stochastic process given by a one-sided distribution Lp[x). Let us consider for the 
moment a Markovian case for which the density of still being one-sided, has finite moments (e.g. an exponential). 
Approximating ^ by its mean, ^ — )• (^), allows us to evaluate the operational time increment. At = ((^)g(x))“^ At. 
The first equation CD can be discretized as Ax = AryAr^/" and the system CD resolves itself to a single equation with 
a multiplicative noise, dx{t) = i/{xy^°‘g{dt), where ^{x) = {{C)g{x))~^. Since the x—dependence is evaluated from 
the other equation, x(t) corresponds to the initial time for each interval and the ltd interpretation is appropriate. 
The above equation incorporates the medium structure but neglects the memory effects which are important if /? < 1. 
Those effects, in turn, can be approximately taken into account by the subordination of the process x(t) to a random 
time t and, after that decoupling of the medium structure and memory, Eq. CD takes the form 

dx{T) = v{xY^°‘g{dT) 

dt{T) = ^(dr); ( 2 ) 


Eq.(l2]) refers to the case /3 < 1 but v{x) contains (^) corresponding to a Markovian process. The first part of Eq.Q 
describes such a process: a Markovian CTRW where the jump length is governed by and the waiting time is 
Poissonian with a position-dependent rate v{x). The derivation of this equation is presented in Appendix A. The 
system similar to ([2]) (but including a potential) was considered in [I^. It was demonstrated in Ref. 221 that Eq.Q 
corresponds to a fractional kinetic equation [l^ containing two fractional operators and a position-dependent diffusion 
coefficient. The anomalous transport predicted by Eq.Q is described by the variance for the Gaussian case and by 
the fractional moments for a < 2 [^ . 

In the present paper, we analyse consequences on the diffusion process which follow from the dependence of memory 
on the position, defined by Eq. CD- In the case of the Levy flights, we construct a process characterised by a finite 
variance taking into account that the system is finite. The finite variance emerges as a result of a variable noise 
intensity near a boundary and which decreases with the distance. This procedure differs from the usual noise truncation 
because the distribution of g remains unaffected. The paper is organised as follows. In Sec.II we discuss the case 
a = 2 exactly solving Eq. © for the power-law g{x). Sec.Ill is devoted to the Levy flights. We discuss a system 
with the multiplicative noise in various interpretations and infer its general properties (Sec.IIIA). In Sec.IIIB, the 
multiplicative noise serves to model a boundary layer; we derive the variance, which then exists, as a function of time 
and all system parameters. 


II. GAUSSIAN CASE 

We consider a stochastic motion inside a medium containing traps and the particle is subjected to a white noise 
the increments of which obey the Gaussian statistics. The substrate nonhomogeneity is restricted to the time- 
characteristics of the system, determined by g(x), and the dynamics is described by Eq.([T]). We assume the function 
g{x) in a power-law form, 

g{x) ^ |xr (0 > -1). (3) 

This form of the diffusion coefficient was assumed to describe, beside a problem of the diffusion on fractals, e.g. a 
turbulent two-particle diffusion [2^ and transport of fast electrons in a hot plasma . It is encountered in geology 
where a power-law distribution of fracture lengths is responsible for transport in a rock ; this self-similar structure 
of a fracture and fault network is characterised by a fractal dimension and determines a pattern of water and steam 
penetration in the rock [^ . In Eq.(|3D, a negative 9 means that intensity of the random time distribution is largest 
near the origin and diminishes with the distance; it rises with the distance for a positive 9. We look for the density 
distribution of the particle position, p{x,t), by solving Eq.Q. The first equation 0 in the Ito interpretation leads 
to the Eokker-Planck equation 

dpo{x,T) _ d‘^[i^{ x)po{x,T)] 
dx"^ 


dr 


( 4 ) 
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which has the solution 


Po{x, t) = A/'t ^\xf exp(—(2 + O^t) 


( 5 ) 


(A/” denotes a normalisation constant), and corresponds to the Markovian process x(t). The random time, given 
by the second equation (ED, is defined by a one-sided, maximally asymmetric stable Levy distribution L^(t), where 
0 < /3 < 1, and the inverse distribution we denote by /i(t, <). The density p{x,t) results from the integration of the 
densities po{x,t) and h(T,t) over the operational time, 

pOO 

p{x,t)= poix,T)h{T,t)dT. (6) 

do 

To evaluate the integral, it is convenient to use the Laplace transform from Eq.([6D taking into account that h{T,u) = 
exp(—Tu^) [ 2 ^. Then a direct integration yields a Laplace transform from the normalised density. 


p{x,u) = 


\x\i+en 

1^0/2 



( 7 ) 


where Ki^(z) is a modified Bessel function, v = 1/{2 + 0) and c = /3 — /3/(4-|- 20) — 1. The density obtained from the 
inversion of the above transform can be expressed in terms of the Fox H-function, 


p{x,t) 


2 (1-f 6»/2)^+2=//5 
I3t r(-:/) 


(2+e)//3-l„2,0 

^1,2 


\x\<a+e)iP 


(0,1) 

{2 + 0YIH 

(c//3- 

vl2, 1/^), (c//3 -b i//2,1/^) _ 


( 8 ) 


details of the derivation are presented in Appendix B. Eq.® seems complicated but an expression for large \x\ is 
simple: it has a stretched-Gaussian form which follows from an asymptotic formula for the H-function [ 2 ^, 


p{x,t) 


(2 + e)( 2 -l 3 ) exp 




(9) 


where A = (2//3 — l)/[/3^/(^“^/(2 -|- 0 )^A^-P))_ xhe stretched-Gaussian shape of the tail is typical for diffusion on 
fractals Q and emerges in the trap models Q. 

We are interested in the moments of p{x, t); all of them are finite and given by a characteristic function which can 
be derived as a series expansion. Eq. © yields 


°° / I 'in fOO 

p{k, t) = Y^ t)dT, 

t'o (2«)! Jo 


where the moments of po{x, r) directly follow from Eq.®, 

_ _/ r , , n-.2n/(2+e) + l ^[(7 + n + 0)/ {2 + 0)] _„/( 2 +e) 

(Z + U) r[-i /(2 + d)] ^ 

Then the integral resolves itself to the moments of h{T,t) that are given by M 

/ 2n/(2+0) \ ^ r[2n/(2 -I- d) -I- 1] 2n0/(2+9) 

^ T[2nl3/{2 + 0) + l] 

and the final expression for the characteristic function reads 

2 + 0 


p(k,t) = - 


r[-l/(2 + 0)]^„ (2n)! 


( 10 ) 


( 11 ) 


( 12 ) 


H + 2n + «)/(2 +'’)! JlinmJl+l] (n) 


Diffusion properties are determined by the variance the time-dependence of which follows from scaling arguments or 
from Eg. dT^ : (x^)(t) = —d‘^p{k = 0,t)/dk^ ^ f'2P/{2+0). formula indicates both the sub- and superdiffusion. 
Gomparison of the above approximate result with a numerical solution of Eq. CD reveals a reasonable agreement; some 
discrepancies emerge for a small /3 and large 0 22|. 

Finally, we will demonstrate that p{x,t), Eq.®, satisfies the fractional equation 


dp{X,t) p,l -/3 92 g 


( 14 ) 
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where qD] ^ is a Riemann-Liouville operator 




1 d 



fit' 


(15) 


and which equation constitutes a non-Markovian generalisation of a Fokker-Planck equation for CTRW with a position- 
dependent waiting-time distribution [s^. First, we integrate Eq. m over time, 


p{x,t) -Poix) = oD^^-^{\x\~^p{x,t)), (16) 

where poix) is the initial condition, and take the Laplace transform, 

ut^Pix, u) - u^"%(x) = ^(|a;|“®p(a;, u)). (17) 

The differentiation produces a differential equation 

x^p" — 9xp' + [0(1 + 9) — u^x^^yp + = 0 (18) 

and the last term vanishes if poix) = (5(a:). Its particular solution has the form 

Pix,u) = ( ^r+ 9/2 

where /(u) is an arbitrary function. Putting /(u) ^ u'^, we obtain Eq.(|7|. 

The power-law form of gix) suggests an interpretation of the memory position-dependence as a fractal trap structure. 
In this picture, gix) means a density of traps, i.e. the number of traps per unit interval, which is self-similar. The 
density of a fractal embedded in the one-dimensional space equals ^ where df stands for the fractal (Hausdorff) 

dimension [s^ . Comparing the above expression with g{x) given by Eq.(j31) allows us to interpret the parameter 9 by 
relating it to the fractal dimension: 9 = df — 1 for 0 G (—1,0]. Then the expression for the variance takes the form. 


{a:2)(t)^t2/3/(i+d/). 


( 20 ) 


The subdiffusion, observed in the homogeneous case, may turn to the enhanced diffusion if dimension of the trap 
structure is small and the memory sufficiently weak (large /3). The case /3 = 1 is special: (^) does not exist and, for 
d/ = I (0 = 0), the problem resolves itself to the ordinary one-dimensional CTRW which is characterised by a weak 
subdiffusion, (x^it) ~ t/lnt Q. 

The present problem differs from the random walk on fractals when the particle performs jumps according to a 
self-similar pattern @ . The non-Markovian generalisation of that model [3| contains a fractional equation which has 
a form different from Eg. (1141) and the predicted motion is always subdiffusive. The density has a finite value at the 
origin in contrast to the our approach: p{x,t), Eq.([S]), diverges at a: = 0 (0 < 0) where the trapping is strong and the 
evolution, as a function of the physical time, proceeds very slowly near the origin. 


III. LEVY FLIGHTS 

Motion of a massive particle subjected to an additive noise which obeys the Levy statistics with long jumps is 
characterised by the infinite variance; this situation is unacceptable for physical reasons. However, the additive noise 
is an idealisation and in realistic systems the stochastic stimulation may depend on a state of the system. For this 
- multiplicative - noise, the variance may be finite and the anomalous diffusion exponent well determined. In the 
next subsection, we discuss general properties of the non-Markovian Langevin dynamics with that noise; the memory 
effects are taken into account by a subordination of the dynamical process to the random time. One can expect that 
the multiplicative noise emerges near a boundary where the environment structure is more complicated than in the 
bulk. The boundary effects modelled in terms of the multiplicative noise will be discussed in Subsection IIIB. 


A. Multiplicative noise 

The generalisation of Eq. m including the multiplicative noise is the following 

dxir) = fix)gidT) 
dtir) = gix)^idT) 


( 21 ) 
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and, in the case of the first equation (I21L we must decide how this equation is to be interpreted, namely at which 
time f{x{t)) is to be evaluated. More precisely, one defines a stochastic integral as a Riemann integral, 

/ f[x{T)]dT]{T) f[{l - Xi)x{Ti-i) + \ix{n)\[r]{Ti) - ri{Ti-i)\, (22) 

“'o i=l 

where the interval (0,t) has been divided in n subintervals (n —>■ oo). The parameter 0 < A/ < 1 determines the 
interpretation and corresponds, in particular, to the ltd (II) (A/ = 0), Stratonovich (SI) (A/ = 1/2) and anti-Ito 
interpretation (All) (A/ = 1). II applies, in particular, if the noise consists of clearly separated pulses, e.g. for a 
continuous description of integer processes and is used in the perturbation theory due to its simplicity. It is well-known 
that for the Gaussian processes the ordinary rules of the calculus are valid in the case of SI, in contrast to the other 
interpretations, which allows us to transform the equation with the multiplicative noise into an equation with the 
additive noise by a simple variable change. As regards the general stable processes, we observe the same property: 
the numerical solutions of the Langevin equation by using Ea. (l22l) with A/ = 1/2 agree with those where a standard 
technique of the variable change is applied [3^ . Obviously, the ordinary rules of the calculus are always valid if one 
defines the white noise as a limit of the coloured noise ; equivalence of that limit and SI is an important property 
of the multiplicative Gaussian processes. In the case of Eq. d^Tl) . we obtain the equation with the additive noise by 
changing the variable, y{x) = ■^jj, in the first of Eq. m- Next, we decouple the medium structure and memory 

(cf. Eq.®) and get the equation dyi^r) = i'{x{y)y^°ri{dT). It contains a multiplicative noise in II and corresponds to 
the Eokker-Planck equation [s^ 


dpoiy,T) ^ d°‘[i^{x{y))po{y,T)] 
dr 5|?/|“ 

Solution in the original variable x directly follows from the solution of Eq. 

Po{x,t) = -^po{y{x),T), 

fix) 


(23) 


(24) 


and the density as a function of the physical time - determined by the equation dt{T) = ^{dr) - results from Eq.®. 
In the following, we restrict our considerations to a power-law form of gix): 


'(x) = 


I — 9 a 




(25) 


similar to Eq.®. Then the waiting time is either large near the origin corresponding to a large intensity of the 
random time density {9 < 0) or the probability of long rests rises with the distance (0 > 0). Eg. (1^51) can be solved if 
one neglects terms higher than |fc|“ in the characteristic function [solls^ which, after inverting the Fourier transform, 
yields tail of the density (IMl) : po{y,T) it corresponds to the Levy-stable asymptotics [s^. The backward 

transformation of the variable, y ^ x, produces the asymptotics po{x,t) f{x)~^y{x)~^~°‘, which indicates that 
finite moments can exist. In particular, the variance is finite if f{x) satisfies the condition 


lim 

X—^OC 


f{x)y{x)^+°‘ 


= 0 . 


(26) 


We assume from now on that /(x) has the power-law form. 


fix) = |a;| (27) 

which allows us to obtain exact solutions. Then Ea. (l26|) implies the condition for the hnite variance: 7 > 2/a — 1. 

Inserting the functions g{x) and fix), given by Eg. (1251) and Eg. (|27)) . to Ea. (l2T]) we obtain after elimination of the 
position-dependent factor in the subordination equation and straightforward calculations the following set of the 
Langevin equations, 

dyir) = D|?/r®/(^+^^ 77 (dr) 

dtir) = ^idr), (28) 

where D = (1 + 7 )“®/(i+'i'). The first of the above equations corresponds to the Eokker-Planck equation, 

dPoiy,x) _ na d°‘[\y\-°‘^/’-^+'^'^poiy,T)] 

Ot ^ ^ 
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which in the limit of small wave numbers is satisfied by a density determined by the following characteristic function 
in respect to the variable y = Y^|a;|^''''’'sign x, 

Poik, r) « 1 - {Aorfo |A:|“; (30) 

in the above equation, eg = 1/(a+ 0/(1 +7)) and 

Eg. (1501) corresponds to the stable density if 0 S (—(1 + 7), (1 + l)/ot) [1^. The final density p{x,t) results from 
the transformation to the physical time by means of Eq.(j6l) which we rewrite as a Fourier transform in respect to y. 
Using Eq. dSOl) yields 


poo 

p{k,t) = 1 — / T^^h{T,t)di 

Jo 


and the integral can be easily evaluated if we express the function h(r, i) by the H-function [3$ 


1 


h{T,t) = 

13t 


T-1//3 

(1,1) 

t 

(l,l//3) . 


(32) 


(33) 


Changing the variable in the integral (155]) . ^ jt, we get 

poo 

p{k, t) = l- (34) 

where H{^) means the H-function in Eq. ([33|l . We evaluate the resulting Mellin transform and obtain the Fourier 
transform corresponding to the stable distribution for small \k\. Inversion of this Fourier transform produces the final 
expression, 


Pix,t) = A-H-^^^\xrL^it-^-^\yix)\/A), 

where we applied Ea. (l24|l and denoted 

A = Alo/‘^r[l + ce]i/“/r[l + /3ce]i/“. 

The solution (1551) has the asymptotics 


p(x,t) 


(35) 


(36) 


(37) 


The variance can be directly evaluated from (|55)l : 

pOO pOO 

{x'^){t) = 2 x^p{xA)dx = 2A-H-P^o x^+^L^{A-\l + -i)-^t-f^^0x^+^)dx. (38) 

Jo Jo 

Change of the variable and evaluating the Mellin transform yields the final expression, 

+ 7)l=«->r (-^) r (i + yT^) (y^) (39) 

which means that motion is always subdiffusive. 

On the other hand, if one understands the noise p in Ea. (l21l) in a sense of II, the Langevin equation for x(t) contains 
a multiplicative term and the counterpart of Ea. (l29l) reads 

dpo{x,T) _ i9°[|a;|~ °^^+T'Vo(cc,T)] 
d\x\°‘ 


dr 


(40) 
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FIG. 1: (Colour online) Distributions obtained from a numerical integration of Ea. (l21ll with g{x) = \x\^°‘ and f{x) = for 
a = 3/2, 9 = 2/3, j3 = 1/2, 7 = 1/2 and t = 0.1. Curves correspond to the following values of A/: 0 (magenta), 0.1 (blue), 0.2 
(cyan), 0.5 (red) and 1 (green) (from bottom to top in the upper part and from top to bottom in the lower part). The straight 
lines correspond to a power-law function with the index 2.5, 2.9 and 3.25 (from top to bottom). 


Solving the above equation in the diffusion limit and evaluation of the integral ([S]) produces the distribution p(x, t) 
with the asymptotics ~ Comparison with Ea. (l37t indicates a qualitative difference between II and SI for 

the Levy flights: whereas for II presence of the multiplicative factor in the noise influences only the time-dependence, 
for SI it modifies the tail shape and makes finiteness of the variance possible. In the face of that difference, it is 
interesting to check what the density distributions for the other interpretations look like. To find those distributions, 
we have to resort to the numerical analysis. According to Ea. (l22l) . the discretized form of the Langevin equation for 
the arbitrary interpretation is given by the following expression, 

— Xji -j- [(1 Xj'jXn ~t“ XjXn^l] ^ 

where rjn is sampled from a symmetric stable distribution according to a well-known algorithm and h is a time 
step. The expression for Xn+i is not explicit and Eq. m can be exactly solved only for a few values of 7 ; in general, 
it must solved numerically at every integration step. For that purpose, we applied the parabolic interpolation scheme 
(the Muller method) [i^. A simple modification of the standard method [ 2 ^ allows us to evaluate the variable x 
as a function of the physical time - determined by the second equation (ED - without an explicit derivation of the 
subordinator ^(t). Examples of distributions for a few values of A/ are presented in Fig.l, separately for the central 
part and for the tails. We observe a similar slope of the tail for all A/ > 0.2; in particular, it is almost identical for 







SI and AIL Near the origin, in turn, some differences emerge and the height of the peak rises with A/. For all the 
interpretations, p{0,t) — 0 which is a consequence of the divergence of f{x) in the origin. 


B. Boundary effects and anomalous diffusion 


Let us assume that the particle is subjected to a noise the intervals of which are governed by the symmetric Levy 
stable distribution. The transport proceeds inside a medium with traps and intensity of the random time density 
is given by the function g{x); the dynamics is described by the Langevin equation with the additive noise, Eq.(IT|). 
Such systems are characterised by the infinite variance which, for massive particles, violates physical principles and 
attempts were undertaken to suppress long tails in the Levy density. A simple remedy is to introduce a modification of 
the stable distribution to make the tail steeper. Such a truncation may be assumed as a simple cut-off [4l| or involve 
some rapidly falling function: e.g. an exponential [i^ or a power-law where f}>2 — a ji^. Processes involving 

the truncated distributions actually converge to the normal distribution, according to the central limit theorem, but 
the power-law tails may be visible for a long time due to a slow convergence. On the other hand, variance becomes 
finite if one takes into account a finiteness of the particle velocity (Levy walk) [l| . In the present approach, the finite 
variance results from a variable intensity of the noise in the region close to the boundary whereas intervals of the 
noise are always distributed according to the stable distribution. As a consequence, the Langevin equation acquires 
a multiplicative noise and the requirement of the finiteness of the variance imposes a condition on the form of the 
position-dependent noise intensity. We assume, in addition, that the boundary effects do not affect the trap structure. 
Presence of that noise is natural: one can expect that the additional complication of the environment structure in 
the vicinity of the boundary introduces a dependence of the noise on the process value and requires a more general 
approach than an equation with the additive noise. Emergence of the multiplicative noise near a boundary has been 
experimentally demonstrated for some physical systems. For example, a description of the colloidal particles diffusion 
in terms of a constant diffusion coefficient appears possible only if a particle remains far from any boundary (d^ . 
Moreover, presence of the multiplicative noise in a description of particles near a wall is necessary to reach a proper 
thermal equilibrium (^ . 

We assume that the boundary effects become important at a distance \x\ = L and then the Langevin equation 
acquires the multiplicative noise. Its intensity is parametrised by f{x) as a falling power-law function, 


J I for |a:| < L 

L'*\x\~'^ for ja:j > L, 


(42) 


where 7 > 0; the dynamics is governed by Eq. m and the noise in the first equation will be interpreted according to 
SI. A new variable. 


y{x) 


X for |a:| < L 

T^[l+ (sign X for ixi > L, 


(43) 


allows us to get rid of the multiplicative factor f{x) in Eq. d^Tl) and the equation corresponding to the first equation 
(051) takes the form 


dx{T) = |a:| ^r]{dT) for |a;| < L 

dyix) = [{\y\ — L){1 +for |a:| > L. (44) 


In general, the density distribution can be obtained in a closed form only for |a:| < L and \x\ ^ L. We consider, 
at the beginning, the case 9 — Q corresponding to a uniform trap distribution and follow a similar method as in the 
preceded subsection: solve the Fokker-Planck equation and integrate over the operational time. The final result for 
\x\ > L reads 


p(x, t) = A-H-PI^L-^\xrL^(t-P/^\y{x)\/A), (45) 

where A = F(I -|- /3)“^/“. The asymptotic form of the above equation, p{x,t) ^ reveals the meaning of 

the parameter 7 : the limit 7 —>■ 00 , for which f{x) = I — 0(|x| — L), yields p{x, t) —>■ 0 (|x| > L), i.e. x = ±L becomes 
an absorbing barrier and the surface is reduced to single points. For a finite value of 7 , the system is not strictly 
confined but the probability density of finding the particle in the outer region rapidly decreases with the distance if 7 
is large. Moreover, for a large L, a very long time is needed to observe the particle at distances markedly larger than 
L. The diffusion problem is well-defined since the finite variance exists if 7 > 2/a — 1 and, in the following, we will 
calculate the variance on this assumption. The system defined by Eg. (1441) changes its properties at \x\ = L and one 
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can expect a different diffusion behaviour in the surface region, compared to the bulk. If L is sufficiently large, the 
dynamics resolves itself to the truncated Levy flights and the diffusion properties are equivalent to the Gaussian case 
[dfil ]. providing the time is relatively small; then one gets the standard anomalous diffusion law, 

(46) 


Variance in the limit t ^ oo follows from a direct evaluation of the integral. One can easily demonstrate that the 
contribution from |a;| < L is small for a large time. Then 




2^-/3/a 

AL~> 





(47) 


and, introducing a new variable x' = t we obtain 


y = L + {t oo) 


1 + 7 

for any x' which allows us to reduce Eq. d+fl) to the form 

2^2/3/a(l+7) 


1+7 


(X^)(t) = 


AL'r 




L-^ 


A(1 +7)‘ 


1 +^ dx, 


(48) 


(49) 


where xq = ^ 0. The integral can be evaluated by applying the standard properties of the H-function 

and, after lengthy but straightforward calculations, we obtain the final formula. 


{x^){t) = L“T'"^(1 + 7 )“"^r(l +^)-"rr(-C 7 )r(ac 7 )sin (50) 

7ra \ 1 + 7 y 

where c-y = 2la(\ + 7 ). Since < 1, the motion is always subdiffusive: the variance in the limit t —>■ 00 rises with 
time slower than linearly and also slower than in the case of a small time, Eg. (1461) . The slope explicitly depends on 
a, in contrast to Ea. (l46l) . and drops to zero for a sharp edge (7 ^ 00 ). Therefore, we observe two diffusion regimes; 
the standard form of the variance (HSl) . which is typical for the Gaussian case and the truncated Levy flights, emerges 
for relatively small times and corresponds to trajectories abiding not far from the origin. On the other hand, when at 
large time the surface region becomes important, diffusion is weaker. Those two diffusion regimes are illustrated in 
Fig.2 for the Gauchy distribution {a = 1) where the variance was obtained bom p{x,t) by a direct integral evaluation: 
slopes agree with Ea. (l46)l and Eo. (1511)1 . 


If d 7 ^ 0, Ea. (l44ll is not manageable analytically except a non-physical case |a;| ^ L. Then the variance was 
determined by a numerical solving of Eg. (!?!]) where the multiplicative noise was treated according to Eg. (HD) with 
A/ = 1/2. The numerical calculations show that the time-dependence of the variance has the same form as Eo. (l46l) 
but with a modified index, ~ This form was found for all the parameters - 0, a and /3 - if L was large. The above 
observation is illustrated in Fig.3 where some examples of {x^)it) are presented. Dependence of the slope on the 
parameters is presented in the subsequent figures. Fig.4 shows /i as a function of 9 for two values of a; it diminishes 
with 9 and the numerical results reveal a dependence p = 0.5/(l + c{a)9) (the coefficient c(a) is indicated in the 
figure). The dependence of /x on a is presented in Fig.5 for two values of 9, both negative and positive. Whereas in 
the former case yt is almost constant - and larger than the value predicted by Ea. (l46|) - for the positive 9 we observe 
a linear growth in a wide range of a; yia) becomes flat only at large a. Finally, Fig .6 shows that the slope for a 
positive (negative) 9 rises with /3 weaker (stronger) than for 0 = 0 and both dependences are linear. The superdiffusion 
emerges when 9 is negative and /3 large. 

The proposed formalism is applicable to diffusion problems where memory is connected with nonhomogeneous 
medium structure. Therefore, we conclude with some remarks concerning the parameters in the above analysis; they 
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FIG. 2: (Colour online) Variance as a function of time for a = 1 and (5 = 1/2. The curves (from bottom to top) correspond 
to: 1.7 = 4 and L = 10; 2. 7 = 2 and L = 10; 3. 7 = 2 and L = 100. Straight-line segments (marked by the red lines) on the 
left-hand side have the slope 1/2 and those on the right-hand side: 1/5, 1/3 and 1/3. 



FIG. 3: (Colour online) Variance as a function of time for a few sets of the parameters a, 9 and /?; the other parameters: 
L = 10 and 7 = 2 . 


may be useful to compare the results with an experiment. We interpreted the multiplicative noise f{x) in Eg. (1211) 
according to SI. This interpretation is distinguished in stochastic problems because it constitutes a white-noise limit 
of coloured noises for any a [s^. However, the other interpretations may also be important. For example, the 
experimental analysis of the colloidal particles diffusion near the boundary favours All (A/ = 1) [i^, as we have already 
mentioned. We performed the numerical calculations for All, similar to those for SI, and found the same diffusion 
properties in respect both to the exponent and the proportionality coefficient. This conclusion could be anticipated 
from Fig.l: tails of the distribution for both interpretations are very similar. Moreover, results presented in Fig.4-6 
are independent of 7 and L if L is sufficiently large. The other parameters have a straightforward interpretation; a 
defines the jump statistics, /3 is responsible for the trapping time characterising the depth of the effective trapping 
potential and 9 is responsible for the trap distribution. 
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FIG. 4: Slope of the time-dependence of the variance, t'', as a function of 6 for = 0.5 and two values of a: 0.5 (lower points) 
and 1.5 (upper points). The other parameters; 7 = 2 and L = 10. Solid lines mark the function /r = 0.5/(l -|- c{a)6), where 
c(0.5) = 1.51 and c(1.5) = 0.51. 



FIG. 5: Slope as a function of a for two values of the other parameters are the same as in Fig.4. Solid line marks the function 
At = 0.2-I-0.12a. 


IV. SUMMARY AND CONCLUSIONS 

The stochastic motion in a medium with traps was studied in terms of the Langevin equation and the anomalous 
diffusion exponent was determined. The memory effects were taken into account by a subordination of a Markovian 
process to a physical, random time and the nonhomogeneity of the trap structure by a dependence of the time lag 
on the position, modelled by a non-negative function g(x). If one decouples effects related to the trap structure 
and memory, the problem resolves itself to a multiplicative process subordinated to the random time. The Langevin 
equation in the operational time describes a Markovian jumping process with a variable rate of the Poissonian 
waiting-time distribution. The density distribution can be exactly derived if g{x) has a power-law form. It has a 
stretched-Gaussian shape and this result is similar to the prediction of the quenched trap model ii- However, in 
contrast to that model, the variance may rise with time faster than linearly; beside the subdiffusion, observed for the 
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FIG. 6: Slope as a function of /3 for two values of 6, the other parameters are the same as in Fig.4. Solid lines mark the 
functions /r = 0.09 + 0.545/1 and /r = —0.05 + 1.276/1 for the positive and negative 6, respectively. 


homogeneous case, the enhanced diffusion emerges. If one interprets the position-dependence of the memory in terms 
of a variable trap density and assume it in a power-law form, such a density corresponds to a fractal pattern. Then 
the fractal dimension determines the diffusion properties: the smaller the dimension, the faster the variance grows 
with time. 

Levy flights are characterised by the infinite variance but a finite size of the system may modify the random 
stimulation. In contrast to the usual truncation procedure where long jumps are eliminated, our approach imposes 
a restriction on the noise intensity by introducing a multiplicative noise to the Langevin equation. This noise may 
follow from a complicated medium structure near the boundary which requires a generalisation of a simple modelling 
in terms of the additive noise. The resulting density distributions have fast falling tails. Therefore, diffusion inside the 
substrate is well-determined and may be quantified in terms of the variance the time-dependence of which does not 
depend of the system size. Moreover, the diffusion properties appear robust in respect to a particular interpretation 
of the multiplicative noise in the Langevin equation. On the other hand, diffusion inside a layer near the boundary is 
weaker than in the bulk and we observe two regimes with a different anomalous diffusion exponent /r. This exponent 
reflects both the magnitude of the memory in the system, described by the parameter j3, and the position-dependence 
of the memory, given by the parameter 0: it diminishes with 6 and rises with /3. In contrast to the case 0 = 0, 
we observe a dependence of /r on a but only for small a and positive 0. This conclusion points at a subtle relation 
between the nonhomogeneous distribution of the random time statistics and the noise statistics in respect to the 
diffusion properties of systems with the Levy flights. 


APPENDIX A 

In this Appendix, we demonstrate that the first part of Eq.(l2]) traces back to a Markovian jumping process [^. 
This stationary process is defined by a jump-size distribution Q{x) in a form of the Levy a-stable and symmetric 
distribution with a characteristic function, 

Q{k) = exp(-A“|A:|“) (a ^ 1, AT > 0). (Al) 

The particle performs instantaneous jumps and then rests for a time given by a waiting-time distribution which is 
Poissonian, 

w{t) = (A2) 

where v{x) denotes a position-dependent rate. The transition probability for infinitesimal time intervals At reads 

Ptrix, /S.t\x', 0) = [1 — v{x')]lS.t5{x — x') + Q{x — x')v{x')At, (A3) 
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where the first term corresponds to the case that no jump occurred within At and the second one that exactly one 
jump occurred. The differentiation over time, 




lim 
—^0 


J Ptr{x, At\x' ,0)p{x' ,t)dx' 


p{x,t) 


/At, 


(A4) 


produces a master equation: 


■^P{x,t) = -v(x)p{x,t) + 


Q(x', x)v{x')p{x', t)dx'. 


(AS) 


In the diffusion limit A: —>■ 0, Eq. (EH can be approximated by Q{k) w 1 — K‘^\k\°‘ and inserting this expression into 
the Fourier transformed Ea. dASI) yields [30l| 

=-K°‘\k\°‘F[v{x)p{x,t)]. (A6) 

Inversion of the above equation yields a fractional Fokker-Planck equation in the form 

dp{x,t) ^^^d°^[v{x)p{x,t)\ 

Sr- = —ajJF—’ ‘ ’ 

where the Weyl-Riesz operator is defined by the inverse Fourier transform: |fc|“). On the other hand, 

Ea. (IA7l) follows from the first of Eq.([21) in the Ito interpretation [ 3 ^ . 


APPENDIX B 


In Appendix B, we derive Eq.(|51). The Bessel function Kjy(z) can be expressed in terms of the Fox H-function [47 
which formula, in our case, reads 


K,(^) 


2 ^ 0.2 


T 


{-v/2,l),{v/2,l) 


(Al) 


After applying standard properties of the H-function and straightforward calculations we obtain the density in the 
form 


p{x,u) 


I {l+e/2Y+'^‘=/P 

^ f(p/7) 


{2+e)/P-lTT2,0 

^ 0,2 




{c/fi-v/2,\/fi),{c/fi-rv/2,\/fY 


(A2) 


where £, = {2 + 9) = ku. Inversion of the transform enhances the order of the H-function [4^; in the 

case of the function in Ea. (IA2l) . the inversion formula yields 


-H 


2,0 

1,2 


( 0 , 1 ) 

(c//3 - v/2, 1//3), (c//3 + u/2, 1//3) 


(A3) 


Finally, we obtain Eq.([5])- 
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